Option pricing with log-stable Levy processes 
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Abstract 

We model the logarithm of the price (log-price) of a financial asset as a random variable 
obtained by projecting an operator stable random vector with a scaling index matrix E_ onto a 
non-random vector. The scaling index E_ models prices of the individual financial assets (stocks, 
mutual funds, etc.). We find the functional form of the characteristic function of real powers of 
the price returns and we compute the expectation value of these real powers and we speculate on 
the utility of these results for statistical inference. Finally we consider a portfolio composed of 
an asset and an option on that asset. We derive the characteristic function of the deviation of the 
portfolio, sj;, defined as a temporal change of the portfolio diminished by the the compound 
interest earned. We derive pseudo-differential equations for the option as a function of the log- 
stock-price and time and we find exact closed-form solutions to that equation. These results 
were not known before. Finally we discuss how our solutions correspond to other approximate 
results known from literature, in particular to the well known Black & Scholes equation. 
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1 Introduction 

Early statistical models of financial markets assumed that asset price returns are independent, identically 
distributed (iid) Gaussian variables. [1]. However, evidence has been found [2] that the returns exhibit 
power law (fat) tails in the high end of the distribution. Except at very high frequencies or short times ([2]), 
a better statistical description for many financial assets is provided by a model where the logarithm of the 
price is a heavy tailed one-dimensional Levy /x-stable process [3, 4, 5, 6]. Since the tail parameter jj, that 
measures the probability of large price jumps will vary from one financial asset to the next, a model based on 
operator stable Levy processes [7] is appropriate. This model allows the tail index to differ for each financial 
asset in the portfolio. Hence we formulate a model where the log-price is a projections of an operator stable 
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random vector onto a predefined direction (this projection determines the portfolio mix). The cumulative 
probability distribution of the log-price diminishes as a mixture of power laws and thus the higher-order 
moments of the distribution may not exist and the characteristic function of the distribution may not be 
analytic. 

Due to the constraints on the size of this paper we only include new results leaving proofs for further 
publications. 

2 The model of the stock market 

In this section we define the model. In the following we recall certain known properties of operator stable 
distributions and we derive Fourier transforms of real powers of projections of operator stable vectors onto 
a non-random vector. In subsections (2.2) and (2.3) we derive Fourier transforms of operator stable random 
vectors for particular forms of parameters of the distribution. 

2.1 The basic properties and new results 

Let log(St) be the logarithm of the price of the portfolio (log-price) at time t. We assume that the temporal 
change of the log-price is composed of two terms, a deterministic term and a fluctuation term viz: 

dS t St+dt — St _ f , . 

— = = adt + a-dL t (1) 

ot Of 

The parameters a G M. (the drift) and the elements of the D dimensional vector a := (<ti, . . . ,<Jd) (the 
portfolio mix) are assumed to be non-random constants. The random vector L t is (strictly) operator stable, 
meaning that it is an operator- normalized limit of a sum of some independent, identically distributed (iid) 
random vectors Xi. We have 

\nt\ 

L t := lim n~ E V X t (2) 

n — >oo A — ' 

i=l 

where n~ E = cxp(— E Inn) and E is a real Z?-dimcnsional matrix such that the equality holds in distribution. 
The class of distributions of the former vectors related to a given matrix E is termed an attraction domain 
of an operator stable law. Members of such class are usually unknown. 

We now recall some known facts [7] concerning the operator stable probability density and its Fourier 
transform ui(k) := Ts[uj](k). 

The following identities hold: 

uj z (x) = lu Ei (t- E x) dct(t- E ) for all t > 0. (3) 

and 

d)*(fc) = exp(-^(fc)) = u{t ET k) = Cj Et {k) (4) 

where E T is the transpose of E and <p is the negative logarithmic characteristic function of the random 
vector L\. In the following we assume that that function is even: 

4>{k) = 4>(-k) (5) 
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The identities (3) are termed as a self-similar property of the random walk L t . 

A motivation for introducing model (1) is statistical inference of parameters of a distribution that de- 
scribes real financial data. In this context it is useful to know analytically the distribution of a real power 
of the integrated fluctuation term in (1). In general this is not known. Here we derive some new results for 
operator stable Levy distributions. Denote by vf\z) the pdf of a random process := (a ■ L^j and by 
vf^ (k) its Fourier transform. For (3 > 1 the identity holds: 



j d\w(t ET a\)K,W («, A) for f3 ^ 1 

-OO 

Cj{t ET ak) for /? = 1 



(6) 



where 



/C (/3) (M) := 



OO 

^ / dte- kt (te- tXtt + r^e 4 ^*) if [/3\ is even 



if L/9J is odd 



(7) 



and a := a/ \a\ and c :— exp(i7r/(2/3)), \i :— exp(— iir{/3} / (f3)). The symbols [/3\ and {/?} mean the biggest 
integer not larger then (3 and the fractional part of [3 respectively, and 6 := ka 13 . 
In addition for even values of [f3\ we have: 



7T 



dze 



sin(Q 



1//3 



r 



where the integration line Z(r, /?)) reads: 

i(t, )3) := [0, r^oo] U [0, (-1)^00] 



(8) 



(9) 



The identities (8) and (9) may be useful for describing the magnitude of the fluctuations of a random walk. 

Thus it follows that the fractional moments of the scalar product are obtained by differentiating the 
Fourier transform at k = 0. We will obtain closed form results for these moments in section (2.4). Here we 
only recall that in the non-Gaussian case, due to (5), we have: 



E 



00 when (3 is even 
when (3 is odd 

In addition for (3 > there the moment exists only if (3 does not exceed a certain threshold value. 



(10) 



It is our objective to price options on the portfolio of stocks driven by operator stable fluctuations (see 
section (3)). By this we mean a theory that 1) allows inference of the the stable index E and the Levy 
measure of the whole vector of stock prices (market)from a statistical sample and 2)hcdges against risk in 
the market by the construction of an appropriate option. To the best of our knowledge, this has not yet 
been achieved. 

The generic properties of operator stable probability distributions and their marginals are described in 
[7]. Here we recall some known facts and we analyse two particular cases of the stable index. The Fourier 
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transform Lu(k) is uniquely determined via the stable index E and the log-characteristic function <fi confined 
to a unit sphere. This can be seen by representing the vector k in the Jurek coordinates viz 

k = rf.e R (11) 
= 1. Using the scaling relation (4) we get: 

exp(-r^(^)) = u r *{6 % ) = u{rf } ■ % ) = Co(k) = exp(-0(fc)) (12) 
and thus 

4>{k) = r % 4>{6 % ) (13) 

Since it follows from the Jordan decomposition theorem that every matrix E is, in a certain basis, a block 
diagonal matrix the set of all possible jump intensities u is narrowed down to few classes of solutions only, 
each one corresponding to a particular Jordan decomposition of the matrix E. We now firstly investigate 
a few classes of solutions as a function of E and subsequently the generic solution for an arbitrary E. The 
existence results in this field are given in [7]. We stress that, contrary to [7], we aim at computing the 
characteristic functions and the fractional moments in closed form rather than only showing their existence. 



2.2 Pure scaling 



In this case E = (D/j.) 1 I where I is a D dimensional identity matrix and n > is a constant. From (11) 
we see that: 

k 



k = rf D ^6 % 



(14) 



hence 



and 6 = kj 



and so 



u)(k) = cxp 



The /3-marginal probability density function from (6) reads: 



d\ exp <^ -t| A| D ^0(— sign(A)) \ /C (/3) {ka , A) 



(15) 



(16) 



where the kernel JC is defined in (7). From the properties of the Gamma function we obtain easily the 
fractional moment of the scalar product as: 



E 



(a ■ L t f]= COS, D M ) • (at^y . ( cos ^ 



e 2 



^/(£>m) 



where 



W,Dy) := 



( 2^ ^(3 + 1 T{-^) 



r(- 



-)■ 



r(-#) 



(17) 



(18) 



Here <f>± := <fr(±). The moment exists for (3 < Dfi. The prcfactor (18) in (17) fits in with the known 
result for the fractional moment of a modulus of a stable variable (see equation (3.6) page 32 in [9]). For the 
derivations of that result by means of the Mellin-Stieljes transform see [25, 26] and by means of characteristic 
functions see [27]. 
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2.3 Scaling & rotation 

In this case D = 2 and we chose: 



£ = 



& (2m)- 1 



(19) 



, cos(/3) — sin(/3) , 

Clearly the trace Tr[_E] = \x . We denote by Op := | | a two dimensional rotation by 

sin(/3) cos(/3) 

an angle /3 e M. The mapping: 



blog(r) 



A; G 



(20) 



changes the length of k by a multiplicative factor A 2 ^ 1 and rotates the vector by an angle — &log(r). The 

2fi 



Jurek coordinates read: = 



and 0£ = O26 ( uiog(|fe|)(fc/ 



and so 



<D(fc) = cxp 



2 M . / fc 
t>2bti log(|fe|) 



The /3-marginal probability density function and the fractional moment read: 



u\ \k) 



dAcxp (-t|A| 2 ^(0 2 ^ log( | A |)^)} JC^ika^ A) 



(21) 



(22) 



and 



27T 

(^it)'] =C0S, 2^)- (vt^y. ( cos (^)e^?) ~ | driP'WiO,*)- (23) 



respectively Here £(/?, 2/i) is defined in (18). The moment exists for < 2)i. Since the moment depends on 
the average of a power of the log-characteristic function over the unit sphere we conclude that the knowledge 
of the moments does not determine the distribution in a unique manner. 

In the following section we compute the fractional moments of the scalar product in the generic case of 
a operator stable distribution. 



2.4 The generic case 

Assume that the stable index has D different eigenvalues {A p }^ =1 that are either real or pairwise complex 
conjugate. Then the following spectral decomposition holds: 

jf =£-Diag(A) -QT 1 (24) 

where Diag(A) := and such that the matrix O is unitary 

£ • £ f = 0? ■ Q, = 1 ( 25 ) 

Then from (24), from the definition of the operator r& := cxp(E r log(r)) and from the Cayley-Hamilton 
theorem we easily arrive at the identity: 

r I T :=£-Diag(r A ) ■ QT 1 (26) 
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where Diag(r A ) := (5ijr Xj )®._ 1 From (11) and (26) we obtain following equations for the Jurck coordinates 
t '■= Ta\ and r := of the vector crA. We denote Xj := Qj + iTj and we have: 

A =(£i§C) =y=i 1 and l = *(p i Qu% > J =^=- ( 27 ) 

where di := 0~ 1 &j are projections of the unit vector a onto the eigenvectors of the stable index (rows of 
the matrix O or columns of the matrix O) . If the unit vector is proportional to the Zth eigenvector then 

o i = 5i,i and from (27) we get r = A®;" 1 and 8 r = \- ir '/ @ 'a. 
The fractional moment reads: 



E 



(a-L t f] = £(P,e T i).(\a\t*>f. 

(cos(|/?)e^) i- / d^* (e"**) 1 E aj,^-^,^ < kT „ (e -^)|- ! (28) 
o 3 * J 

Here I is such a number that 6; — Qj > for all j 7^ 2 such that |<jj| > 0, and J := {j \@j = 6;}, D* = 
D — card(J) + 1, and £(/?, 6 ; _1 ) is defined in (18). The moment exists for /39; < 1, meaning if (3 does not 
exceed the inverse biggest real part of eigenvalues of the tail index E_ in the maximal eigenspace containing 
a. We note the following: 

[1] The moment is a product of a real prefactor, the /3th power of the length of the vector cr, a time 
factor r 1 and a complex prefactor. The former prefactor is the same as in the one dimensional case 
whereas the later prefactor is a complex number equal to the support of the random variable (a ■ L^j . In 
particular for symmetric distributions the later factor is real and the support of random variable is given by 

l + (-l)' 3 = cos(|/3)e J f' 3 . 

[1] The integrand in (28) is related to the complement of the J-space, meaning a linear span of eigenvectors 
whose real parts of eigenvalues equal O; . 

[2] The result (28) is in accordance with an existence result (Theorem 8.3.10 in [7]). However we have for 
the first time computed the moment in closed form which will be useful for statistical inference for example 
or for other theoretical work. 

[3] If the multiplicity of 9; is equal to D then card( J) = D and the complement of J is empty and the last 
term in (28) reduces to 

</ e ') := (27T)- 1 j d^ e < (e-*&) (29) 


because the Heaviside function in the integrand is identically equal unity. This is like in the scaling & 
rotation case. 

[4] If <f>{e~^ l &) < 1 then the left hand side of the equality in the subscript of the Heaviside function is positive 
and the Heaviside function may not be identically equal unity. 
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3 The option price 



An option on a financial asset is an agreement settled at time t to purchase (call) or to sell (put) the asset 
at some maturity time T in the future. Here we consider European style options that can be exercised only 
at maturity. This means that boundary conditions are imposed on the option price at maturity t = T. Ex- 
tending the analysis to American style options that may be exercised at any time can be done by considering 
European style options with a different number of exercise times [10] and allowing the number of exercise 
times to go to infinity. 

In order to minimize the risk we now divide the money available between Ns stocks S t and Nc options 
C(S t ; t). The value of the portfolio is then: 



We may, without loss of generality, chose Nc = 1 • 

The portfolio is a stochastic process that is required to grow exponentially with time in terms of its 
expectation value. The rate of growth r is the so-called "riskless" rate of interest and is assumed to be 
independent of time t. 

3.1 Local temporal growth 

Consider the distribution of deviations 



between the interest (e rdt — l)V(t) = (rdt + 0{dt 2 ))V{t) that is earned by the portfolio and the change 
V(t + dt) — V(t) of the price of the portfolio. Does a self-financing strategy exists? Is it possible to choose 
C = C(St;t) subject to a condition Ct = max(SV — K, 0), for some strike price K, such that the expectation 
value of the deviations of the portfolio conditioned on the price of the stock at time t equals zero? Thus we 
require that the deviations have no drift: 



In our model we assume that the above condition is satisfied only for an infinitesimal time change dt 
and is conditioned on the value of the stock price at time t (local temporal growth sec. 3.1). Due to limited 
space we are not able to include a model extension that assumes that the above condition is satisfied for a 
finite dt We will present it in a future publication. 

Our approach is more general than that used in financial mathematics [11] where considerations are 
based on the lack of arbitrage, meaning the assumption that riskless opportunities for making money in 
financial transactions do not exist. We waive that unrealistic assumption and instead require the portfolio 
to increase exponentially with time. 

From equation (1) we have: 



V(t)=N s S t + N c C(St;t) 



(30) 




V(t + dt) - e rdt V(t) 



(31) 



E D [ i \St -0 



(32) 




(33) 
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where we have used the fact that a Levy process is homogeneous in time, meaning that 

L t +dt - L t = L dt (34) 

where = in (34) means an equality is in distribution. We note that (33) is merely a transformation of 
equation (1) and not a solution to that equation. As such equation (33) holds for infinitesimal times dt only. 

From (33) we see that the expectation value of the right hand side conditioned on St and is infinite, 
unless the fluctuations are Gaussian. 

Therefore we will modify the log-characteristic function cf){\) in order to ensure the finiteness of all 
moments. We define: 



&(A) := 0(A) exp(- 



) (35) 



where e > and replace (j) e by (f>. From now on we will work with a fictious process related to the modified 
log-characteristic function, we will solve the option pricing problem for it and at the end of the derivation 
we will take the limit e — > 0. After finishing the derivation we will check analytically if the result ensures a 
risk free portfolio. Firstly we check that conditional expectation value of the stock price is finite. We have: 

E [St +d t - S t \ S t ] = S t (cxp (adt) E [e^'^*)] - l) = S t (exp (adt) e-***-*) - l) < oo (36) 
where in the second equality in (36) we used the following identity: 

E [e(* £ *)] = J dze z 5 (z - $ v 9 . u% (£R = J dze* ■ i- J dke* k ^ ■ v 3 . ldt (£R (37) 

R R R 

= Jdk5(k + t)D ff . £ Jk)= J dk5{k + i)i> 3 . idt {k) = e- dt ^ (38) 

R -i+R 

In the first equality in (37) we inserted a delta function into the definition of the expectation value, in the sec- 
ond equality we used the integral representation of the delta function, in the first equality in (38) we integrated 
over z and £ and we used the integral representation of the delta function in the second equality we shifted the 
integration line by using the Cauchy theorem applied to a rectangle [-R, R] U R + i[0, 1] U — i + [R, R] U R — z[0, 1] 
in the limit R — > oo and in the last equality we used (6) and (4). We make three remarks. Firstly the delta 
function has been analytically continued to complex arguments, ie we have defined it as follows: 

5{k + iq) := jr £^5<p) (jfe) (39) 

Secondly we note that the result (38) holds only for e > because otherwise, all p > 2 terms in the sum (39) 
produce infinite values when integrated with the second term in the integrand. Thirdly we reiterate that it 
is the fictious, modified stock price, related to e > 0, that has a finite expectation value whereas the real 
stock price has of course an infinite expectation value. The option prices that we compute correspond to a 
fictious e-world where stock prices' probabilities have been modified like in (35) and the real option price is 
obtained as a limit of the former as e tends to zero and the sequence of the fictious worlds towards our real 
world. In other words the option pricing problem has indeed no solution in our real world however it has a 
solution in the "complement" of our world by the limit of the e-worlds. 
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We will therefore construct a zero-expectation value stochastic process (31) as a linear combination (30) 
of two stochastic processes St and C(S t ;t) that have both non-zero expectations values. For this purpose 
we will analyze the probability distribution of the deviation variable and work out conditions for the 



option price such that the conditional expectation value E 
deviation of the portfolio: 



is equal zero. Now we compute the 



® (dt) = (V(t + dt) - V{t)) + V(t) (1 - e rdt ) (40) 
= N S (S t+dt - S t ) + (C t+dt - C t ) + V(t) (1 - e rdt ) (41) 
= N s S t (e adt+ff - E «* - l) + (C(S t+dt ;t + dt) - C t ) + V(t) (l - e rdt ) (42) 

= N s S t [e adt+3 - L ^ - 1 ) + ^dt + — { (e adt+3 - Ldt - l) +V(t)(l- e rdt ) (43) 

m— 1 

In (43) we have expanded the price of the option in a Taylor series to the first order in time and to all 
orders in the price of the option. In that we have assumed that the price of the option is a perfectly smooth 
function of the price of the stock. This may limit the class of solutions. In particular, solutions may exist, 
where the price of the stock is a function satisfying the Holder condition: 

\C(S t+d t;t + dt)-C t \ < A\S t+d t-S t \ A (44) 

for any St+dt an d S t , a constant A and a Holder exponent A G [0, 1) and thus the price of the option can be 
expanded in a fractional Taylor series [14] in powers of S t + d t — S t . We will seek for these solutions in future 
work. 

The process is a sum of infinitely many terms that have non-zero expectation values. We could 

compute its expectation value directly using (36) and re-sum the series. However, we will instead calculate 
the characteristic function of the process conditioned on the value of the process St at time t. This 

means that we propagate the process St by an infinitesimal value dt and we compute the characteristic 
function of the increment and we require the zero value derivative of the characteristic function to be equal 
zero. This technique is not new, see discussion about solving master equations of Markov processes in 
[17], and it works because of the time-homogeneity of the process (34) and because of the fact that the 
parameters a and a are constant as a function of the process S t . The time- homogeneity follows from the 
infinite divisibility of the process and thus the technique applied here also works in the generic setting of 
Levy processes. We will extend the model according to these lines in future investigations. 

We note that D (dt) in (43) is a function of the scalar product a ■ L dt only and thus the distribution of 
J)^*' is unique functional of the distribution of the scalar product. 

We derive the distribution of ®[ dt ^ now. We define 

D t := (d t Cdt + V(t) (l - e rdt )) (45) 
we condition on the value of the fluctuation a ■ L dt , we use (43), and we get: 

1 + Yl ^ %n{q) cxp (-^-^ - i a )) ( 46 ) 

m=l ' q=0 / 
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where <f> is the negative logarithmic characteristic function of the random vector L\ (see (4)). The expectation 
value of the portfolio deviation conditioned on the value of the price of the stock S t reads: 



where 



E 



® { t dt) \s t 



d{ik) 



k=0 



8,C - rV(.) - (* + -) St (6,, - „) - £ •.gg^ 



dt + O (dt 2 ) 



6 «^ : =E| | (-l) n -m-tqa) and <E r , 
9=0 \ q 



fe=max(n, 2) 



(47) 
(48) 

(49) 



and the log-characteristic function 4>(—iqa) in (48) has been analytically continued to imaginary arguments. 
Here the coefficients ai read: 



n— i -. 

a (*) :=(_!)*-»(*._!)! ^ J] - = (-l) fe "" E 

l<ji<...<j„_i<fe-l g=l ^ 



l<ji<...<j fe _ n <fe-l 9=1 



with = (— l) fe x (fc — 1)!. In addition the coefficients <£„ satisfy: 



(50) 



X>„=0 



(51) 



what follows readily from the fact that J2n=i c4 = for fc > 2. 

In the Gaussian case the coefficients read & n ,ct> = ( — < J ) 2 (n$n,i + n(n — l)fl n 2) and thus (48) yields a 
second order PDE. Since the Levy distribution has been truncated as in (35) and due to (5) the result in 
(49) is real. Indeed the log-characteristic function can be expanded in a Taylor series in even powers of the 
argument only and thus its value at the negative imaginary unit is real. If we did not truncate we would have 
obtained a unrealistic complex result as seen from (15). We reiterate that the limit of truncation threshold 
going to zero (e — > in (35)) will be taken at the end of the calculation only rather than at intermediary 
stages. If we did so at this stage we would have obtained a paradoxical result; an infinite sum of numbers 
<£ n each of which is infinite equals zero. 

The requirement E Df^ \S t — 0(dt) implies a following generalized Black & Scholes equation: 



n=2 



pfn.fi a/~i 

'-ro(^(s)r =rm + {Ns+ os }S ' {6 ^- a) (52) 

In order that we get further insight into the problem, in particular in order that we are able to solve 
equation (52) analytically we find a new expression for the coefficients of the PDE stated in the following 
propositions. 

Proposition 1 The coefficients &k,<t> m (52) read: 



oo 

©M = y^(£)(-l + e-«) fe 



(53) 
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for m € N. Here </>(£) is the inverse Laplace transform of the log-characteristic function of L\ or the Levy 
measure of the process L t . We have: 

4>(0 ■■= ^ / dze^4>(-iza) , <j>(-iz) = J d^e^'fa) (54) 

iR R + 

In the pure scaling case for D = 1 the inverse Laplace transform </> reads: 

where 1^1 ^ = T>± x and the later operators are Marchaud whole axis fractional derivatives One has to bear in 
mind that since the function cp(—iza) may be in general unbounded as a function of z and thus the quantity 
<^>(£) is in general not a function but a functional. 
Proposition 2 The coefficients <E„ in (52) read: 



(-l)/d^(0(e-*-l + £) ifn=l 

R+ (56) 
(-1)" / d^(0£r ifn>l 



From (56) and (55) we see that the coefficients are infinite if D/j, < 2. 

We proceed as follows to solve the PDE (52) . In the definition (56) of the coefficients <B n we truncate 
the upper limit of integration at some threshold value then we solve the generalized Black& Scholes equation 
(52) analytically by Fourier transforming with respect to log(S') and at the end we take the limit of the 
truncation threshold to infinity. Note that this step is essential. Indeed, as seen from (54) and from (35) 
it is not clear if the inverse Laplace transform <fr related to the truncated Levy distribution diminishes fast 
enough away from the origin and thus if the integral in (56) exists. We accomplish this task in section (3.2). 
Prior to doing that we describe how we will compute the number stocks as follows. 

We define a utility function U of the portfolio as a functional of the price of the stock viz: 

T T 

U:= J V(OdZ - J (N s St + C(S; £)) d£ (57) 
o o 

and require (57) to be minimal. We do not investigate here the mathematical subtleties concerned with the 
existence of the stochastic integral (57). The necessary condition is that the variation SU with respect to 
the price of the stock functional is zero. We have: 

T 

SU := J (n 3 + ^§^) SSdt = (58) 
o 

what yields that 

as in the Gaussian case. We note that this choice of the number of stocks ensures the self-financing property 
of the portfolio. Indeed in the Cox-Ross-Rubinstcin binary tree model in discrete time one considers a 
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portfolio composed of a stock and a bond and one derives the number of stocks by requiring contingent 
claim replication, meaning an equality of the portfolio and the claim with probability one (see [28] for 
example). The later result is essentially the same as that in (59). 

Comments We have derived a PDE for the option price that ensures that the derivative of the expectation 
value of the portfolio with compounded interest is zero 

r E[V(t + dt)-e rdt V(t)\S t ] 

hm^o — j t ~ = (60) 

without making any assumptions about the relationship between the drift of the stock price a and the riskless 
rate of interest r. We differ in that from standard models in financial mathematics [20, 21], models that 
assume at the outset that a = r. 



3.2 Final result 



We solve the generalized Black& Scholes equation analytically. Inserting (59) into the second equality in 
(52) we get: 

dC 



d t C = rC-r 



dx 



n=l 



(61) 



where x = \og(S t ). The coefficients (£„ are defined in (56) with the upper limit of integration being truncated 
at some threshold value. Since the coefficients do not depend on x the PDE (61) is converted into a Ordinary 
Differential Equation (ODE) by taking a Fourier transform of the option price with respect to x. This gives: 



where 



and 



d t C{k-t) = (r + H(k))C(k;t) 
C(x;t) := (27T)- 1 J dxC{k- 1 t)e~ tkx 



H{k) :- 



rik + ^2 £n(-ik) r 



rik + V(k) 



(62) 
(63) 

(64) 



We insert (56) into (64) and obtain the following expression for the function V(k) that we call 'the Hamil- 
tonian' after Hagen Kleinert [16]. We have: 

(65) 



V(k) = / d^(£) [(-1 + e-*)(*fc) + e 4fe « - 1] = (-1) 2 / d£/W (0 [(*fce-*) - fcV fc «] 



We see that the integrals in (65) exist. Therefore the limit of the truncation threshold in these integrals 
going to infinity can be performed at this stage. This is what we do now and assume hereafter the whole 
positive real axis in the integration in (65). From (54) we obtain the Hamiltonian 

V(k) = (ik) <j>(-ia) + 4>{-ka) (66) 

Now we come back to equation (62) which we solve subject to an initial condition at maturity as follows: 

C(k; t) = C(k; T) exp {- (r + H(k)) t} (67) 
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where C(k; T) is the Fourier transform of the option payoff C(x; T) at maturity T and r := T — t is the time 
to maturity. This payoff reads: 

{max(e x — K, 0) for a call 
(68) 
maxfif — e^, 0) for a put 

where K is the strike price. The Fourier transform of the payoff is easily computed and it reads: 



for a call 



C(k; T) := f d£C(£; T)e lki = I } ^ *" ; '{ ((,<)) 

1 K + (+Tlki- Tk- 2lTS ( k - 1 )) for a put 

We insert (69) into (67) and invert the Fourier transform for a call. We have: 

C(x;t) f dkC(k;t)e' tkx (70) 

2tt J 

R 

— log(m)+rr — log(m)+rr 

= S t J d^-%r Lr ^ + <t>{-ia)T)-Ke-^ f d & 3 . Zr {ti + <t>{-i5)T) (71) 

— oo — oo 

= S t N l {d l )-Ke- rT N 2 {d l ) (72) 

We recall that here v g j- (^) is the probability density function of the fluctuation term a ■ L T . 
In (72) we have changed the integration variables and simplified the result. Here we defined: 

d d 

Ni(d) := J dee- £ e T *(-^i/ ff .£ T (0 , N 2 (d) := J ^.£ T (0 and di = - log(m) + r (r + 0(-tff)) (73) 

— oo — oo 

In the limit D/i — > 2_ the density z/- ^ (^) in (71) goes into a Gaussian with mean zero and variance 2<tt 
and (72) goes into the Gaussian Black& Scholes equation, see e.g. eqs. (1.6), (1.7) on page 8 in [11]. We end 
this section by stating the price of the portfolio. We have: 

V(t) = N s St + C{x-t) = -Ke- rT N 2 (d 1 ) (74) 

Since the last factor on the right-hand side in (74) depends implicitly on St the unconditional expectation 
value of the portfolio does not increase exponentially as required. Therefore the solution (72) is only an 
approximation. However since the Gaussian Black & Scholes equation is a particular case of (72) it turns 
out that it is also only an approximation. 

The factors in (73) are complex which is of course unrealistic. The reason for that is the following. In 
our approach we assumed that the time change dt is infinitesimally small rather asumming it to be finite 
at the outset and taking the limit dt — > at the end of the calculation. We have checked that the later 
procedure leads to a real result which has the same form as in (72) except that the Levy density Vg.£ (£) goes 
into an inverse Fourier transform of exp (l/2a IJ 'T(tk + fjb)/T{ik)) evaluated at £ — rr<j){—ia) which essentially 
amounts to replacing the expression — k^ + ik by T(ik + n)/T(ik) in some intermediate calculations. The 
inverse Fourier transform in question is essentially equal to the log(SV) — ar process probability density 
function evaluated at the argument £ — rr(j>(—ia). Therefore the price of the option is a discounted present 
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value of the maturity payoff under a risk- neutral probability measure where the measure in question is related 
to the compensated log-price process log(St) — at. Thus we have proven that the risk-neutral option pricing 
method holds in the generic setting of operator stable processes. 

Expressions (73) are difficult to deal with in numerical calculations. Indeed the "typical width" of the 
Fourier transform of the Levy density is cr^r. Since this quantity is small, meaning of the order of 10 -2 for 
stock daily data and for times to maturity of the order of hundreds of days, the use of "primitive" methods 
like Romberg quadratures for evaluating the Fourier integrals requires a very high precision of calculation 
that is much bigger than the precision of the estimated parameters. Therefore we propose to use a more 
sophisticated method for the numerical integrations. This method is described in the Appendix. 

4 Conclusions 

We have applied the technique of characteristic functions to the problem of pricing an option on a stock that 
is driven by operator stable fluctuations. We have developed a technique to ensure that the expectation value 
of the portfolio grows exponentially with time. In doing this we have not, unlike other authors, made any 
assumptions about the analytic properties of the log-characteristic function of the stock price process. Instead 
we have expressed all results in terms of the characteristic function of the operator stable fluctuation L\ . 

Subsequent to successful numerical tests, we ought then to be able to price analytically not only European 
options but also exotic options with a finite number of different exercise times. This should also allow us to 
price American style options by allowing the number of exercise times to become infinite. 

We may also compute the 99th percentile of the probability distribution of the deviation of the portfolio 
(Value at Risk) as a function of a and of the log-characteristic function of the random vector L\ . The 
Value at Risk will be expressed as an integral equation involving the conditional characteristic function of 
the portfolio deviation (46). The resulting integrals will be carried out by means of the Cauchy complex 
integration theorem. 

The results of these calculations will be reported in a future publication. 
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6 Appendix 

We explain how the integrals from the Levy density in (72) are computed numerically in the case D = 1. 
Note that in this case equation (72) can be written as follows: 

C(x; t) = S t N { X (d; z) Ke^N^ (d; z) (75) 

where z := z r + %z% = 4>{—ia)r is a complex number, d := — log(m) + rr is a real number and 

d 

N { ;l(d;z):= J d^ V<yL ^ + z) (76) 

— oo 

for s = 0,1. We note that the factors (76) have a following integral representation that lends itself to 
numerical computations in a straightforward manner. We have: 

d d 

N[ s) L {d;z)= [ d^e-^-— [ dke- lk ^ +z '> ■ e^^l = — [ dke- lkz ■ [ d^e- (s+lk ^ e - T ^ k) (77) 

J 2-7T J A^oo 2-7T J J 

-oo R R —A 

1 r / p — iBd i pt6A \ i f / p — iBd i „t8A\ 

t i- ± -^ L± -)-°-* t) = ^h I *e-.(^^).e-*l7 8 ) 

R R+»s 

j- / d9e- l(e+ls ^ ■ ( 6 ^ j • e -r<K°^s)) (79) 



A^oo 



A^oo 



16 



r<fi(tscr) 



dOe 



sH9( d+Zr )) M i + r m d + Z r)) M2 ■ 

C7 C7 



(80) 



In the first equality in (77) we expressed the Levy stable density through its Fourier transform and in the 
second equality in (77) we changed to order of integration. In the first equality in (78) we integrated over 
£ and we defined — k — is and in the second equality in (78) we shifted the integration line by is in the 
complex plane. In doing this we used the Cauchy theorem applied to a rectangular contour composed of an 
interval [—R,R], of that interval shifted by is and of sections perpendicular to the real axis that complete 
the contour. In the limit R — > oo the integrals over the later sections vanish. In (79) we factorize-d the 
integral and in (80) we performed the limit A — > oo by using the identity: 

o i0A 



liniA- 
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= n8(6) for e 



where 



Mi(6) := J2 e- T ^ p0+ ^ = and M 2 {9) := J2 pe -^(^+ 
p=±i p=±i 

In the pure scaling case in one dimension, from (16), we have 4>(k) — <j)±k 11 and thus: 



»«)) 



(81) 
(82) 



MM := Y, 



-(f>±r(al) ti cos(/j.0 p ) 



cos(a p ) - ism{a p )) and M 2 {6) := J2 V^^^ cos( ^ p) (cos(a p ) - tsin(a p )) 



p=±i 



p=±i 



where 



/. := + s 2 , cos(0 p ) = — ,sin(0 p ) = -, and a p = </>±t(ct/) m sin(/x^ p ) 



(83) 



(84) 



Since, as seen from (82), M X {Q) = 2e- r ^<") and M 2 (0) = the integral in (80) is clearly finite the result 
can be used for numerical calculations. In the Gaussian case /i = 2we have a p = 2(f>±T<r 2 sp9 and thus 



M\{6) := 2cos(20 ±S( 7 2 0T)e^ ±rCT2 ( e2 - s2 ) and M 2 (9) := -2t sin^sa 2 ^)^™ 2 ^ 2 ^ 2 ) 
and z = T(fi(—ia) = z r + iz,- t = —T(f)±a 2 . Inserting (85) into (80) gives: 



(85) 



9 2 2i 



(86) 



= ^ + ^ y (»?) = y d? 7^£ T w = y 



c??7 



2tt 



(87) 



where e = — log(m) + r (r + 0±er 2 (2s — 1)). From (87) we see that the factors coincide with those in the 
Gaussian Black& Scholcs formula. 
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